home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Programming / Source / IP_Pvc / Pvc_Object.m < prev    next >
Encoding:
Text File  |  1992-05-10  |  21.9 KB  |  1,045 lines
  1. #import <sys/message.h>
  2. #import <mach_error.h>
  3. #import <soundkit/Sound.h>
  4. #import <string.h>
  5. #import "Pvc_Object.h"
  6.  
  7. void bitreverse( float *buf,int N );
  8.  
  9. @implementation Pvc_Object : Object
  10. {
  11.  
  12. /*
  13.  
  14.   int        N,
  15.         N2,
  16.         Nw = 0,
  17.         Nw2, 
  18.         decim = 0,
  19.         interp = 0,
  20.         in,
  21.         on,
  22.         valid = -1,
  23.         eof = 0;
  24.   float     freqmlt = 0.,
  25.         srate = 44100.,
  26.         pi,
  27.         twopi,
  28.         synt = 0.,
  29.         *Hwin,
  30.         *Wanal,
  31.         *Wsyn,
  32.         *input,
  33.         *buf,
  34.         *channel,
  35.         *output;
  36.   BOOL        obank = 0,
  37.         aflag = 0,
  38.         sflag = 0;
  39.   NXZone    *ozone;
  40.  
  41. */
  42.  
  43. }
  44.  
  45. + create {
  46.     id         newInstance;
  47.     newInstance = [ self new ];
  48.     return newInstance;
  49. }
  50.  
  51. - init {
  52.     [super init];
  53.     [self setDefaults];
  54.     inputFilename = malloc(1024);
  55.     status = IDLE;
  56.     return self;
  57. }
  58.  
  59. - setInputFilename: (char *) arg
  60. {
  61.     int err;
  62.     if (inputFilename && !(strcmp(arg,inputFilename)))
  63.         return self;
  64.     if (insfh != NULL)
  65.     SNDFree(insfh);        
  66.     if (strlen(arg) > 1023) 
  67.         return nil;
  68.     strcpy(inputFilename,arg);
  69.     err = SNDReadSoundfile(inputFilename, &insfh);
  70.     if (err != SND_ERR_NONE)
  71.     return nil;
  72.     return self;
  73. }
  74.  
  75. - setOutputFilename: (char *) arg
  76. {
  77.     outputFilename = arg;
  78.     return self;
  79. }
  80.  
  81. static BOOL isPowerOfTwo(int n)
  82.     /* Query whether n is a pure power of 2 */
  83. {
  84.     while (n > 1) {
  85.     if (n % 2) break;
  86.     n >>= 1;
  87.     }
  88.     return (n == 1);
  89. }
  90.  
  91. - checkArgs:(char *)msg 
  92.  /* If failure, sets msg to error message and returns nil.
  93.     Msg should be at least 1024 long */
  94. {
  95.  
  96.     if (Nw == 0)
  97.     Nw = N;
  98.  
  99.     if (interp == 0)
  100.       interp = decim;
  101.  
  102.     obank = freqmlt != 0.;
  103.     N2 = N>>1;
  104.     Nw2 = Nw>>1;
  105.     if (interp > Nw) {
  106.     sprintf(msg,"I must be less than or equal to window size.");
  107.     return nil;
  108.     }
  109.     if (decim > Nw) {
  110.     sprintf(msg,"D must be less than or equal to window size.");
  111.     return nil;
  112.     }
  113.     if (decim > Nw) {
  114.     sprintf(msg,"D must be less than or equal to window size.");
  115.     return nil;
  116.     }
  117.     if (!isPowerOfTwo(N)) {
  118.     sprintf(msg,"FFT size (N) must be a power of 2.");
  119.     }
  120.     return self;
  121. }
  122.  
  123. - allocateDataspace {
  124.  
  125.     int        width,
  126.         size;
  127.  
  128.     if ( ozone != NULL )
  129.       NXDestroyZone(ozone);
  130.  
  131.     ozone = NXCreateZone( 10 * Nw * sizeof(float), vm_page_size, 0 );
  132.  
  133. /* analysis window */
  134.     Wanal = (float *) zspace( ozone, Nw, sizeof(float) );
  135.  
  136.  /* synthesis window */
  137.     Wsyn = (float *) zspace( ozone, Nw, sizeof(float) );
  138.  
  139. /* input buffer */
  140.     input = (float *) zspace( ozone, Nw, sizeof(float) );
  141.  
  142. /* hamming window */
  143.     Hwin = (float *) zspace( ozone, Nw, sizeof(float) );
  144.  
  145. /* FFT buffer */
  146.     buf = (float *) zspace( ozone, N, sizeof(float) );
  147.  
  148. /* analysis channels */
  149.     channel = (float *) zspace( ozone, N+2, sizeof(float) );
  150.  
  151. /* output buffer */
  152.     output = (float *) zspace( ozone, Nw, sizeof(float) );
  153.  
  154.     SNDGetDataPointer( insfh, (char **)&inData, &insize, &width );
  155.     if (outsfh != NULL)
  156.         SNDFree(outsfh);
  157.     SNDAlloc( &outsfh, 
  158.          ( (int) ((N+insize) * width * ( (float) interp / decim ))),
  159.             SND_FORMAT_LINEAR_16, insfh->samplingRate,
  160.             insfh->channelCount, 0);
  161.     SNDGetDataPointer(outsfh, (char **)&outData, &size, &width );
  162.  
  163.     srate = insfh->samplingRate;
  164.  
  165.     return self;
  166. }
  167.  
  168. - (SNDSoundStruct *)outsfh
  169. {
  170.   return outsfh;
  171. }  
  172.  
  173. - (SNDSoundStruct *)insfh
  174. {
  175.   return insfh;
  176. }  
  177.  
  178. - setDefaults {
  179.  
  180.   N = 1024;
  181.   Nw = 0;
  182.   valid = -1;
  183.   decim = 128; 
  184.   interp = 0;
  185.   eof = 0;
  186.   freqmlt = 0.;
  187.   srate = 44100.;
  188.   pi = 4 * atan(1.);
  189.   twopi = 8 * atan(1.);
  190.   synt = 0.;
  191.   obank = 0;
  192.   inPoint = 0;
  193.   outPoint = 0;
  194.  
  195.   return self;
  196. }
  197.  
  198. - (float)howFar {
  199.   return ((float)outPoint) / ((outsfh->dataSize)/sizeof(short));
  200. }
  201.  
  202. - N: (int) aN {
  203.     N = aN;
  204.     return self;
  205. }
  206.  
  207. - Nw: (int) aNw {
  208.     Nw = aNw;
  209.     return self;
  210. }
  211.  
  212. - decim: (int) aDecim {
  213.     decim = aDecim;
  214.     return self;
  215. }
  216.  
  217. - interp: (int) aInterp {
  218.     interp = aInterp;
  219.     return self;
  220. }
  221.  
  222. - freqmlt: (float) aFreqmlt {
  223.     freqmlt = aFreqmlt;
  224.     return self;
  225. }
  226.  
  227. - srate: (float) aSrate {
  228.     srate = aSrate;
  229.     return self;
  230. }
  231.  
  232. - synt: (float) aSynt {
  233.     synt = aSynt;
  234.     return self;
  235. }
  236.  
  237. - aflag: (BOOL) aAflag {
  238.     aflag = aAflag;
  239.     return self;
  240. }
  241.  
  242. - sflag: (BOOL) aSflag {
  243.     sflag = aSflag;
  244.     return self;
  245. }
  246.  
  247. - writeOutputSound {
  248.  
  249.   int err;
  250.   err = SNDWriteSoundfile( outputFilename, outsfh );
  251.   if (err != SND_ERR_NONE)
  252.       return nil;
  253.   return self;
  254. }
  255.  
  256. - kill
  257. {
  258.     [self resume];
  259.     eof = 1;
  260.     return self;
  261. }
  262.  
  263. - resume
  264.   /* Gets called from the other thread */
  265. {
  266.     kern_return_t ec;
  267.     msg_header_t msg =    {0,                   /* msg_unused */
  268.                            TRUE,                /* msg_simple */
  269.                            sizeof(msg_header_t),/* msg_size */
  270.                            MSG_TYPE_NORMAL,     /* msg_type */
  271.                            0};                  /* Fills in remaining fields */
  272.     if (eof != 2 && appToObjPort != PORT_NULL)
  273.       return nil;
  274.     msg.msg_local_port = PORT_NULL;
  275.     msg.msg_remote_port = appToObjPort;
  276.     ec = msg_send(&msg, SEND_TIMEOUT, 0);
  277.     if (ec == SEND_TIMED_OUT)
  278.       ; /* message queue is full, don't need to send another */
  279.     return self;
  280. }
  281.  
  282. - pause
  283. {
  284.     eof = 2;
  285.     return self;
  286. }
  287.  
  288. - (int)status
  289. {
  290.     return status;
  291. }
  292.  
  293. - _waitForMessage
  294.     /* This is the main loop for a separate-threaded performance. */
  295. {
  296.     struct {
  297.         msg_header_t header;
  298.         char data[MSG_SIZE_MAX];
  299.     } msg;
  300.     status = PAUSED;
  301.     (void)port_allocate(task_self(), &appToObjPort);
  302.     msg.header.msg_size = MSG_SIZE_MAX;
  303.     msg.header.msg_local_port = appToObjPort;
  304.     (void)msg_receive(&msg.header, MSG_OPTION_NONE,0);
  305.     (void)port_deallocate(task_self(),appToObjPort);
  306.     appToObjPort = PORT_NULL;
  307.     status = RUNNING;
  308.     return self;
  309.   }
  310.  
  311. /* phase vocoder */
  312.  
  313. - runPvc
  314.     /* Returns self if successfully writes file, else nil. */
  315. {
  316.  
  317.     status = RUNNING;
  318. /* initialize input and output time values (in samples) */
  319.     in = -Nw;
  320.     if ( decim )
  321.     on = (in * interp) / decim;
  322.     else
  323.     on = in;
  324.  
  325. /* main loop--perform phase vocoder analysis-resynthesis */
  326.  
  327.     while ( eof != 1) {
  328.  
  329.         if (eof == 2)  {
  330.       [self _waitForMessage];
  331.       if (eof == 1)
  332.         break;
  333.     }
  334.         
  335. /* increment times */
  336.     in += decim;
  337.     on += interp;
  338.  
  339. /* analysis: input D samples; window, fold and rotate input
  340.    samples into FFT buf; take FFT; and convert to
  341.    amplitude-frequency (phase vocoder) form */
  342.  
  343.     eof = [self shiftin];
  344.     [self fold];
  345.     [self rfft:YES];
  346.     [self convert];
  347.  
  348. /* at this point channel[2*i] contains amplitude data and
  349.    channel[2*i+1] contains frequency data (in Hz) for phase
  350.    vocoder channels i = 0, 1, ... N/2; the center frequency
  351.    associated with each channel is i*f, where f is the
  352.    fundamental frequency of analysis R/N; any desired spectral
  353.    modifications can be made at this point: pitch modifications
  354.    are generally well suited to oscillator bank resynthesis,
  355.    while time modifications are generally well (and more
  356.    efficiently) suited to overlap-add resynthesis */
  357.  
  358. /* oscillator bank resynthesis */
  359.  
  360.     if ( obank ) {
  361.  
  362.         [self oscbank]; 
  363. /* oscbank( channel, N2, srate, Nw, interp, freqmlt, output ); */
  364.         [self shiftout];
  365.     } 
  366.  
  367. /* overlap-add resynthesis */
  368.  
  369.     else {
  370.         [self unconvert];
  371.         [self rfft:NO];
  372.         [self overlapadd];
  373.         [self shiftout];
  374.     }
  375.     }
  376.     if (![self writeOutputSound]) {
  377.     status = IDLE;
  378.     return nil;
  379.     }
  380.     status = IDLE;
  381.     return self;
  382. }
  383.  
  384.  
  385. /* overlap-add FFT analysis/resynthesis without phase unwrapping */
  386.  
  387. - runFFT
  388. {
  389.  
  390. /* argument checking */
  391.     status = RUNNING;
  392.  
  393. /* initialize input and output time values (in samples) */
  394.     in = -Nw;
  395.     if ( decim )
  396.     on = (in * interp) / decim;
  397.     else
  398.     on = in;
  399.  
  400. /* main loop--perform phase vocoder analysis-resynthesis */
  401.  
  402.     while ( !eof ) {
  403.  
  404. /* increment times */
  405.     in += decim;
  406.     on += interp;
  407.  
  408. /* analysis: input D samples; window, fold and rotate input
  409.    samples into FFT buf; take FFT; and convert to
  410.    amplitude-frequency (phase vocoder) form */
  411.  
  412.     eof = [self shiftin];
  413.     [self fold];
  414.     [self rfft:YES];
  415.     [self convertFFT];
  416.     
  417. /* at this point channel[2*i] contains amplitude data and
  418.    channel[2*i+1] contains frequency data (in Hz) for phase
  419.    vocoder channels i = 0, 1, ... N/2; the center frequency
  420.    associated with each channel is i*f, where f is the
  421.    fundamental frequency of analysis R/N; any desired spectral
  422.    modifications can be made at this point: pitch modifications
  423.    are generally well suited to oscillator bank resynthesis,
  424.    while time modifications are generally well (and more
  425.    efficiently) suited to overlap-add resynthesis */
  426.  
  427. /* overlap-add resynthesis */
  428.  
  429.     [self unconvertFFT];
  430.     [self rfft:NO];
  431.     [self overlapadd];
  432.     [self shiftout];
  433.       }
  434.     status = IDLE;
  435.     return self;
  436. }
  437.  
  438. /* shift next D samples into righthand end of array A of
  439.    length N, padding with zeros after last sample (A is
  440.    assumed to be initially 0); return 0 when more input
  441.    remains, otherwise return 1 after N-2*D zeros have been
  442.    padded onto end of input */
  443.  
  444. - (int) shiftin
  445. {
  446.   register int     i;
  447.  
  448.   if ( valid < 0 )        /* first time only */
  449.     valid = Nw;
  450.  
  451.   for ( i = 0 ; i < Nw - decim ; i++ )
  452.     *(input+i) = *(input + i + decim);
  453.  
  454.   if ( valid == Nw ) {
  455.     for ( i = (Nw - decim); i < Nw; i++) {
  456.       if ( inPoint >= insize ) {
  457.           *(input+i) = (float) (*(inData + inPoint) / 32768.);
  458.     valid = i - (Nw-decim);
  459.     break;
  460.       }
  461.       *(input+i) = (float) (*(inData + inPoint) / 32768.);
  462.       ++inPoint;
  463.     }
  464.   }
  465.   if ( valid < Nw ) {        /* pad with zeros after EOF */
  466.     for (i=valid; i < Nw; i++)
  467.       *(input+i) = 0.;
  468.     valid -= decim;
  469.   }
  470.   return (valid <= 0);
  471. }
  472.  
  473.  
  474. /* if output time n >= 0, output first I samples in
  475.    array A of length N, then shift A left by I samples,
  476.    padding with zeros after last sample */
  477.  
  478. - shiftout
  479. {
  480.   int    i;
  481.  
  482.   if ( on >= 0 ) {
  483.     for ( i=0; i < interp; i++ )
  484.       *(outData + outPoint + i) = (short) (32767. * *(output+i));
  485.       
  486.     outPoint += interp;
  487.   }
  488.  
  489.   for ( i = 0; i < Nw - interp; i++ )
  490.     *(output+i) = *(output + i + interp);
  491.  
  492.   for ( i = (Nw - interp); i < Nw; i++ )
  493.     *(output+i) = 0.;
  494.  
  495.   return self;
  496. }
  497.  
  498. /* S is a spectrum in rfft format, i.e., it contains N real values
  499.    arranged as real followed by imaginary values, except for first
  500.    two values, which are real parts of 0 and Nyquist frequencies;
  501.    convert first changes these into N/2+1 PAIRS of magnitude and
  502.    phase values to be stored in output array C; the phases are then
  503.    unwrapped and successive phase differences are used to compute
  504.    estimates of the instantaneous frequencies for each phase vocoder
  505.    analysis channel; decimation rate D and sampling rate R are used
  506.    to render these frequency values directly in Hz. */
  507.  
  508. - convert
  509. {
  510.   static int     first = 1;
  511.   static float     *lastphase,
  512.         fundamental,
  513.         factor;
  514.   float     phase,
  515.         phasediff;
  516.   int         real,
  517.         imag,
  518.         amp,
  519.         freq;
  520.   float     a,
  521.         b;
  522.   int         i;
  523.  
  524. /* first pass: allocate zeroed space for previous phase
  525.    values for each channel and compute constants */
  526.  
  527.     if ( first ) {
  528.       first = 0;
  529.       lastphase = (float *) space( N2+1, sizeof(float) );
  530.       bzero(lastphase,(N2+1)*sizeof(float));
  531.       fundamental = srate / (float) (N2<<1);
  532.       factor = srate / (decim * twopi);
  533.     } 
  534.  
  535. /* unravel rfft-format spectrum: note that N2+1 pairs of
  536.    values are produced */
  537.  
  538.     for ( i = 0; i <= N2; i++ ) {
  539.       imag = freq = ( real = amp = i<<1 ) + 1;
  540.       a = ( i == N2 ? *(buf+1) : *(buf+real) );
  541.       b = ( i == 0 || i == N2 ? 0. : *(buf+imag) );
  542.  
  543. /* compute magnitude value from real and imaginary parts */
  544.  
  545.       *(channel+amp) = hypot( a, b );
  546.  
  547. /* compute phase value from real and imaginary parts and take
  548.    difference between this and previous value for each channel */
  549.  
  550.       if ( *(channel+amp) == 0. )
  551.     phasediff = 0.;
  552.       else {
  553.     phasediff = ( phase = -atan2( b, a ) ) - *(lastphase+i);
  554.     *(lastphase+i) = phase;
  555.     
  556. /* unwrap phase differences */
  557.  
  558.     while ( phasediff > pi )
  559.       phasediff -= twopi;
  560.     while ( phasediff < -pi )
  561.       phasediff += twopi;
  562.       }
  563.  
  564. /* convert each phase difference to Hz */
  565.  
  566.       *(channel+freq) = (phasediff * factor) + (i * fundamental);
  567.     }
  568.     return self;
  569. }
  570.  
  571. - convertFFT
  572. {
  573.  
  574.   int         amp,
  575.         phase;
  576.   float     a,b;
  577.   register int    i;
  578.  
  579.  
  580.     for ( i = 0; i <= N2; i++ ) {
  581.       phase = (amp = i<<1) + 1;
  582.       a = ( i == N2 ? *(buf+1) : *(buf+amp) );
  583.       b = ( i == 0 || i == N2 ? 0. : *(buf+phase) );
  584.  
  585. /* compute magnitude value from real and imaginary parts */
  586.  
  587.       *(channel+amp) = hypot( a, b );
  588.       *(channel+phase) = -atan2( b, a );
  589.     }
  590.     return self;
  591. }
  592.  
  593.  
  594. /* unconvert essentially undoes what convert does, i.e., it
  595.   turns N2+1 PAIRS of amplitude and frequency values in
  596.   C into N2 PAIR of complex spectrum data (in rfft format)
  597.   in output array S; sampling rate R and interpolation factor
  598.   I are used to recompute phase values from frequencies */
  599.  
  600. - unconvert
  601. {
  602.   static int     first = 1;
  603.   static float     *lastphase,
  604.         fundamental,
  605.         factor;
  606.   int         i,
  607.         real,
  608.         imag,
  609.         amp,
  610.         freq;
  611.   float     phase;
  612.  
  613. /* first pass: allocate memory and compute constants */
  614.  
  615.     if ( first ) {
  616.     first = 0;
  617.     lastphase = (float *) space( N2+1, sizeof(float) );
  618.     fundamental = srate / (float) (N2<<1);
  619.         factor = twopi * interp / srate;
  620.     } 
  621.  
  622. /* subtract out frequencies associated with each channel,
  623.    compute phases in terms of radians per I samples, and
  624.    convert to complex form */
  625.  
  626.     for ( i = 0; i <= N2; i++ ) {
  627.     imag = freq = ( real = amp = i<<1 ) + 1;
  628.  
  629.     if ( i == N2 )
  630.         real = 1;
  631.  
  632.     *(lastphase+i) += *(channel+freq) - i * fundamental;
  633.     phase = *(lastphase+i) * factor;
  634.     *(buf+real) = *(channel+amp) * cos( phase );
  635.  
  636.     if ( i != N2 )
  637.         *(buf+imag) = - *(channel+amp) * sin( phase );
  638.     }
  639.     return self;
  640. }
  641.  
  642.  
  643. - unconvertFFT
  644. {
  645.  
  646.   int         amp,
  647.         phase;
  648.   register int    i;
  649.  
  650.  
  651.     for ( i = 0; i <= N2; i++ ) {
  652.     phase = (amp = i<<1) + 1;
  653.  
  654.     if ( i == N2 )
  655.       *(buf+1) = *(channel+amp) * cos( *(channel+phase) );
  656.     else
  657.       *(buf+amp) = *(channel+amp) * cos( *(channel+phase) );
  658.  
  659.     if ( i != N2 )
  660.         *(buf+phase) = - *(channel+amp) * sin( *(channel+phase) );
  661.     }
  662.     return self;
  663. }
  664.  
  665.  
  666. /* make balanced pair of analysis (A) and synthesis (S) windows;
  667.    window lengths are Nw, FFT length is N, synthesis interpolation
  668.    factor is I, and osc is true (1) if oscillator bank resynthesis 
  669.    is specified */
  670.  
  671. - makewindows
  672. {
  673.   int i;
  674.   float sum;
  675.  
  676. /* basic Hamming windows */
  677.  
  678.   for ( i = 0; i < Nw; i++ )
  679.     *(Hwin+i) = *(Wanal+i) = *(Wsyn+i) = 0.54 - (0.46 * cos( 
  680.         twopi * i / (Nw - 1) ));
  681.  
  682. /* when Nw > N, also apply interpolating (sinc) windows to
  683.    ensure that window are 0 at increments of N (the FFT length)
  684.    away from the center of the analysis window and of I away
  685.    from the center of the synthesis window */
  686.  
  687.     if ( Nw > N ) {
  688.      float x;
  689.  
  690. /* take care to create symmetrical windows */
  691.  
  692.     x = -(Nw - 1) / 2.;
  693.     for ( i = 0; i < Nw; i++, x += 1. )
  694.         if ( x != 0. ) {
  695.         *(Wanal+i) *= N * sin( (pi * x) / N ) / (pi*x);
  696.         if ( interp )
  697.            *(Wsyn+i) *= interp * sin( (pi*x) / interp ) / (pi*x);
  698.         }
  699.     }
  700.  
  701. /* normalize windows for unity gain across unmodified
  702.    analysis-synthesis procedure */
  703.  
  704.     for ( sum = i = 0; i < Nw; i++ )
  705.     sum += *(Wanal+i);
  706.  
  707.     for ( i = 0; i < Nw; i++ ) {
  708.      float afac = 2./sum;
  709.      float sfac = Nw > N ? 1./afac : afac;
  710.     *(Wanal+i) *= afac;
  711.     *(Wsyn+i) *= sfac;
  712.     }
  713.  
  714.     if ( Nw <= N && interp ) {
  715.     for ( sum = i = 0; i < Nw; i += interp )
  716.         sum += *(Wsyn+i) * *(Wsyn+i);
  717.     for ( sum = 1./sum, i = 0; i < Nw; i++ )
  718.         *(Wsyn+i) *= sum;
  719.     }
  720.     return self;
  721. }
  722.  
  723.  
  724. /* input I is a folded spectrum of length N; output O and
  725.    synthesis window W are of length Nw--overlap-add windowed,
  726.    unrotated, unfolded input data into output O */
  727.  
  728. - overlapadd
  729. {
  730.   int     i,
  731.     n;
  732.  
  733.   n = on;
  734.   while ( n < 0 )
  735.     n += N;
  736.   n %= N;
  737.  
  738.   for ( i = 0 ; i < Nw ; i++ ) {
  739.     *(output+i) += *(buf+n) * *(Wsyn+i);
  740.     if ( ++n == N )
  741.       n = 0;
  742.   }
  743.   return self;
  744. }
  745.  
  746.  
  747. /* multiply current input I by window W (both of length Nw);
  748.    using modulus arithmetic, fold and rotate windowed input
  749.    into output array O of (FFT) length N according to current
  750.    input time n */
  751.  
  752. - fold
  753. {
  754.  
  755.   int     i,
  756.     n;
  757.  
  758.   n = in;
  759.  
  760.   for ( i = 0; i < N; i++ )
  761.     *(buf+i) = 0.;
  762.  
  763.     while ( n < 0 )
  764.     n += N;
  765.  
  766.     n %= N;
  767.     for ( i = 0; i < Nw; i++ ) {
  768.       *(buf+n) += *(input+i) * *(Wanal+i);
  769.       if ( ++n == N )
  770.     n = 0;
  771.     }
  772.     return self;
  773. }
  774.  
  775. /* If forward is true, rfft replaces 2*N real data points in buf with
  776.    N complex values representing the positive frequency half of their
  777.    Fourier spectrum, with *(buf+1) replaced with the real part of the Nyquist
  778.    frequency value.  If forward is false, rfft expects buf to contain a
  779.    positive frequency spectrum arranged as before, and replaces it with
  780.    2*N real values.  N MUST be a power of 2. */
  781.  
  782. - rfft: (BOOL) forward
  783. {
  784.   float     c1,c2,
  785.           h1r,h1i,
  786.         h2r,h2i,
  787.         wr,wi,
  788.         wpr,wpi,
  789.           temp,
  790.         theta;
  791.   float     br,bi;
  792.   int         i,
  793.         i1,i2,i3,i4,
  794.         N2p1;
  795.  
  796.     theta = pi/N2;
  797.     wr = 1.;
  798.     wi = 0.;
  799.     c1 = 0.5;
  800.     if ( forward ) {
  801.     c2 = -0.5;
  802.     [self cfft:forward];
  803.     br = *buf;
  804.     bi = *(buf+1);
  805.     }
  806.     else {
  807.     c2 = 0.5;
  808.     theta = -theta;
  809.     br = *(buf+1);
  810.     bi = 0.;
  811.     *(buf+1) = 0.;
  812.     }
  813.     wpr = -2. * pow( sin( 0.5*theta ), 2. );
  814.     wpi = sin(theta);
  815.     N2p1 = (N2 <<1) + 1;
  816.     for ( i = 0; i <= N2 >> 1; i++ ) {
  817.     i1 = i<<1;
  818.     i2 = i1 + 1;
  819.     i3 = N2p1 - i2;
  820.     i4 = i3 + 1;
  821.     if ( i == 0 ) {
  822.         h1r =  c1 * ( *(buf+i1) + br );
  823.         h1i =  c1 * ( *(buf+i2) - bi );
  824.         h2r = -c2 * ( *(buf+i2) + bi );
  825.         h2i =  c2 * ( *(buf+i1) - br );
  826.         *(buf+i1) =  h1r + (wr * h2r) - (wi * h2i);
  827.         *(buf+i2) =  h1i + (wr * h2i) + (wi * h2r);
  828.         br =  h1r - (wr * h2r) + (wi * h2i);
  829.         bi = -h1i + (wr * h2i) + (wi * h2r);
  830.     }
  831.     else {
  832.         h1r =  c1 * ( *(buf+i1) + *(buf+i3) );
  833.         h1i =  c1 * ( *(buf+i2) - *(buf+i4) );
  834.         h2r = -c2 * ( *(buf+i2) + *(buf+i4) );
  835.         h2i =  c2 * ( *(buf+i1) - *(buf+i3) );
  836.         *(buf+i1) =  h1r + wr*h2r - wi*h2i;
  837.         *(buf+i2) =  h1i + wr*h2i + wi*h2r;
  838.         *(buf+i3) =  h1r - wr*h2r + wi*h2i;
  839.         *(buf+i4) = -h1i + wr*h2i + wi*h2r;
  840.     }
  841.     wr = ((temp = wr) * wpr) - (wi * wpi) + wr;
  842.     wi = (wi * wpr) + (temp * wpi) + wi;
  843.     }
  844.     if ( forward )
  845.     *(buf+1) = br;
  846.     else
  847.     [self cfft:forward];
  848.  
  849.     return self;
  850. }
  851.  
  852. /* cfft replaces float array x containing NC complex values
  853.    (2*NC float values alternating real, imagininary, etc.)
  854.    by its Fourier transform if forward is true, or by its
  855.    inverse Fourier transform if forward is false, using a
  856.    recursive Fast Fourier transform method due to Danielson
  857.    and Lanczos.  NC MUST be a power of 2. */
  858.  
  859. - cfft: (BOOL) forward
  860. {
  861.   float     wr,wi,
  862.         wpr,wpi,
  863.         theta,
  864.         scale;
  865.   int         mmax,
  866.         ND,
  867.         m,
  868.         i,j,
  869.         delta;
  870.  
  871.   ND = N2<<1;
  872.   bitreverse( buf, ND );
  873.  
  874.   for ( mmax = 2; mmax < ND; mmax = delta ) {
  875.     delta = mmax<<1;
  876.     theta = twopi / ( forward? mmax : -mmax );
  877.     wpr = -2. * pow( sin( 0.5*theta ), 2. );
  878.     wpi = sin(theta);
  879.     wr = 1.;
  880.     wi = 0.;
  881.     
  882.     for ( m = 0; m < mmax; m += 2 ) {
  883.       register float rtemp, itemp;
  884.       
  885.       for ( i = m; i < ND; i += delta ) {
  886.     j = i + mmax;
  887.     rtemp = ( wr * *(buf+j) ) - ( wi * *(buf+j+1) );
  888.     itemp = ( wr * *(buf+j+1) ) + ( wi * *(buf+j) );
  889.     *(buf+j) = *(buf+i) - rtemp;
  890.     *(buf+j+1) = *(buf+i+1) - itemp;
  891.     *(buf+i) += rtemp;
  892.     *(buf+i+1) += itemp;
  893.       }
  894.       wr = ((rtemp = wr) * wpr) - (wi * wpi) + wr;
  895.       wi = (wi * wpr) + (rtemp * wpi) + wi;
  896.     }
  897.   }
  898.  
  899. /* scale output */
  900.  
  901.   scale = forward ? 1./ND : 2.;
  902.   
  903.   { 
  904.     register float *bi=buf, *be=buf + ND;
  905.     
  906.     while ( bi < be )
  907.       *bi++ *= scale;
  908.   }
  909.   return self;
  910. }
  911.  
  912. /* bitreverse places float array x containing N/2 complex values
  913.    into bit-reversed order */
  914.  
  915. void bitreverse( float *buf,int N )
  916. {
  917.   float     rtemp,itemp;
  918.   int         i,j,
  919.         m;
  920.  
  921.   for ( i = j = 0; i < N; i += 2, j += m ) {
  922.  
  923.     if ( j > i ) {
  924.  
  925.       rtemp = *(buf+j);    /* complex exchange */
  926.       itemp = *(buf+j+1);
  927.       
  928.       *(buf+j) = *(buf+i);
  929.       *(buf+j+1) = *(buf+i+1);
  930.       
  931.       *(buf+i) = rtemp;
  932.       *(buf+i+1) = itemp;
  933.     }
  934.     
  935.     for ( m = N>>1; m >= 2 && j >= m; m >>= 1 )
  936.       j -= m;
  937.   }
  938. }
  939.  
  940.  
  941.  
  942. /* oscillator bank resynthesizer for phase vocoder analyzer
  943.    uses sum of N+1 cosinusoidal table lookup oscillators to 
  944.    compute I (interpolation factor) samples of output O
  945.    from N+1 amplitude and frequency value-pairs in C;
  946.    frequencies are scaled by P */
  947.  
  948. - oscbank
  949. {
  950.   static int     NP,
  951.         tablesize = 8192,
  952.         first = 1;
  953.   static float     Iinv,
  954.         *lastamp,
  955.         *lastfreq,
  956.         *index,
  957.         *table,
  958.            Pinc,
  959.         ffac;
  960.   int         amp,
  961.         freq,
  962.         n,
  963.         chan;
  964.  
  965. /* first pass: allocate memory to hold previous values
  966.    of amplitude and frequency for each channel, the table
  967.    index for each oscillator, and the table itself; also
  968.    compute constants */
  969.  
  970.     if ( first ) {
  971.   
  972.       float twopioL = twopi / tablesize, tabscale;
  973.  
  974.     first = 0;
  975.     lastamp = (float *) space( N2+1, sizeof(float) );
  976.     lastfreq = (float *) space( N2+1, sizeof(float) );
  977.     index = (float *) space( N2+1, sizeof(float) );
  978.     table = (float *) space( tablesize+1, sizeof(float) );
  979.  
  980.     tabscale = ( Nw >= N2 ? N2 : 8 * N2 );
  981.  
  982.     for ( n = 0; n < tablesize+1; n++ )
  983.         *(table+n) = tabscale*cos( twopioL*n );
  984.  
  985.     Iinv = 1. / interp;
  986.     Pinc = freqmlt * ( (float) tablesize ) / srate;
  987.     ffac = freqmlt * pi / ((float) N2);
  988.     if ( freqmlt > 1. )
  989.         NP = (int) ((float) N2 / freqmlt);
  990.     else
  991.         NP = N2;
  992.     }
  993.  
  994. /* for each channel, compute I samples using linear
  995.    interpolation on the amplitude and frequency
  996.    control values */
  997.  
  998.     for ( chan=0; chan < NP; chan++ ) {
  999.  
  1000.       register float     a,
  1001.             ainc,
  1002.             f,
  1003.             finc,
  1004.             address;
  1005.  
  1006.       freq = ( amp = ( chan << 1 ) ) + 1;
  1007.       if ( *(channel+amp) < synt ) /* skip the little ones */
  1008.     continue;
  1009.       *(channel+freq) *= Pinc;
  1010.       finc = ( *(channel+freq) - ( f = *(lastfreq+chan) ) ) *Iinv;
  1011.       
  1012.       ainc = ( *(channel+amp) - ( a = *(lastamp+chan) ) ) *Iinv;
  1013.       address = *(index+chan);
  1014.       
  1015. /* accumulate the I samples from each oscillator into
  1016.    output array O (initially assumed to be zero);
  1017.    f is frequency in Hz scaled by oscillator increment
  1018.    factor and pitch (Pinc); a is amplitude; */
  1019.  
  1020.       for ( n=0; n < interp; n++ ) {
  1021.     *(output+n) += a * *(table + ( (int) address ));
  1022.     address += f;
  1023.     
  1024.     while ( address >= (float) tablesize )
  1025.       address -= (float) tablesize;
  1026.     
  1027.     while ( address < 0. )
  1028.       address += (float) tablesize;
  1029.     
  1030.     a += ainc;
  1031.     f += finc;
  1032.       } 
  1033.  
  1034. /* save current values for next iteration */
  1035.  
  1036.       *(lastfreq+chan) = *(channel+freq);
  1037.       *(lastamp+chan) = *(channel+amp);
  1038.       *(index+chan) = address;
  1039.     }
  1040.     return self;
  1041. }
  1042.  
  1043.  
  1044. @end
  1045.